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ABSTRACT 

In this paper we investigate how convective instabiHties influence heat conduction in the intra- 
cluster medium (ICM) of cool-core galaxy clusters. The ICM is a high-beta, weakly collisional 
plasma in which the transport of momentum and heat is aligned with the magnetic field. The 
anisotropy of heat conduction, in particular, gives rise to instabilities that can access energy 
stored in a temperature gradient of either sign. We focus on the heat-flux buoyancy-driven 
instability (HBI), which feeds on the outwardly increasing temperature profile of cluster cool 
cores. Our aim is to elucidate how the global structure of a cluster impacts on the growth 
and morphology of the linear HBI modes when in the presence of Braginskii viscosity, and 
ultimately on the ability of the HBI to thermally insulate cores. We employ an idealised quasi- 
global model, the plane-parallel atmosphere, which captures the essential physics - e.g. the 
global radial profile of the cluster - while letting the problem remain analytically tractable. 
Our main result is that the dominant HBI modes are localised to the innermost (<20%) re- 
gions of cool cores. It is then probable that, in the nonlinear regime, appreciable field-line 
insulation will be similarly localised. Thus, while radio-mode feedback appears necessary 
in the central few tens of kpc, heat conduction may be capable of offsetting radiative losses 
throughout most of a cool core over a significant fraction of the Hubble time. Finally, our 
linear solutions provide a convenient numerical test for the nonlinear codes that simulate the 
saturation of such convective instabilities in the presence of anisotropic transport. 

Key words: conduction - instabilities - magnetic fields - MHD - plasmas - galaxies; clus- 
ters: intracluster medium. 



1 INTRODUCTION 

Whereas their copious X-ray emission suggests a central cooling 
time much less than the Hubble time, galaxy cluster cores do not 
exhibit cooling flows commensurate with their radiative losses. Ab- 
sent are the mass deposition rates (A/ ^ 10^-10^ Mq yr^^) and 
copious iron line emission that would accompany such flows. In- 
stead, the central temperatures of cool core clusters are typically 
only ~l/3 of the bulk cluster temperatures, while spectroscopi- 
cal ly determined mass depos ition rates are <Q.1M (for a review, 
see lPeterson & Fabian l2006l) . Understanding how these conditions 
are maintained in the face of rapid radiative cooling poses a chal- 
lenge for theorists, who have consequently invoked a plethora of di- 
verse physics to resolve the problem. These include: active galactic 
nucleus feedback, conductive heat transport, and convective turbu- 
lence, whether it be driven by cosmic rays, mergers, or magnetohy- 
drodynamic (MHD) instabilities. 

This paper concerns the contribution of conductive heat trans- 
port to the solution of the cooling flow problem. However, we 
do not take the usual energetics standpoint (i.e. does conduction 
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provide an ample source of heat to offset radiative losses?) but 
rather focus on how conduction generates, and is subsequently 
modified by, plasma instabilities and the MHD turbulence they 
instigate. Because the intracluster medium (ICM) is magnetized 
and weakly collisional it exhibits surprising stabili ty properties 
( lBalbusll200(]l 1200 ll : ISchekochihin et alj|2005l I2OI0I) . In particu- 
lar, the outwardly increasing temperature profile of a cl uster core 
gives rise to a heat-flux buoyancy-driven instability (HBI:^ OuataerJ 
[2OO8) that can drive disordered flows. Numerical simulations show 
that these flows lead to near-complete field-lin e insulation of the 
core, and heat conduction is greatly impeded (iParrish & Ouataerd 
l2008l : iBogdanovic et al.|[2009l : |Parrish, Ouataert fe'sha rma 200^ 
HBI motions not only seem incapable of maintaining the ob- 
serve d temperature profile (as suggested by iBalbus & Revnoldsl 
|2008|) . but actually exacerbate the cooling flow problem. As 
a result, theorists have appealed to turbulent stirring by other 
mean s to reopen the field lines and reinstate the conductive heat- 
flux iP arrish, Quatae rt & Sharmal l20ld : iRuszkowski & OhI l20ld : 
iMcCourt et al . 2010). 

This would bring the issue to a close if it were not for the 
fact that a magnetized, weakly collisional plasma can exhibit ap- 
preciable pressure anisotropy - an effect neglected in most previ- 
ous work. In fact, HBI modes inevitably generate such anisotropies, 



©2011 RAS 



2 Latter & Kunz 



manifested as viscous stresses, and these adversely affect the in- 
stability mechanism. In particular, the wavelengths of the fastest 
growing modes Ahbi increase from very small scales t o those of or - 
der the thermal-pressure scaleheight of the cluster H jKunzll201ll) . 
So, even though the HBI mechanism can be understood with lo- 
cal physics, the instability manifests primarily on global scales and 
may be sensitive to the large-scale structure of the cluster, espe- 
cially via the kinematic viscosity 's strong depe ndence on tempera- 
ture and density. For this reason, iKunj ( 1201 ll) conjectured, on ac- 
count of a WKBJ analysis, that the fastest growing HBI modes 
would favour only the innermost regions of the cluster core, as it 
is here that viscous damping is minimised. If the nonlinear satura- 
tion of the HBI is similarly confined, then the associated field-line 
insulation of the core may be significantly attenuated. 

In this paper, we investigate how pressure anisotropy and the 
global structure of th e ICM influe nce the HBI by generalising the 
linear calculation of iKun j ( 1201 ih . We calculate the HBI modes 
in a simple quasi-global model of the cluster: a plane-parallel at- 
mosphere. The model captures the radial density and temperature 
structure of the ICM, but neglects its angular structure. Though ide- 
alized, this model permits us to test in a transparent way the new 
physics at work, and how the character of the HBI changes in a 
global setting. In addition, it provides a numerical testbed for non- 
linear codes that include anisotropic conduction and viscosity. Our 
main result is that the dominant HBI modes are typically localised 
to the innermost (<20%) regions of cool-core clusters. It is then 
likely that, in the nonlinear regime, appreciable field-line insula- 
tion will similarly localise and become inefficient. These theoret- 
ical results have bee n corrobora t ed by the Braginskii-MHD simu- 
lations presented in iKunz et al.l ( |2012[) . which use domains suffi- 
ciently large to capture this physics. 

The paper is organized as follows. In Section|2]we present the 
model and its governing equations. In Section [3] we derive quasi- 
global equilibrium solutions to these equations, which are then sub- 
jected to linear perturbations in Section |4] There, the HBI eigen- 
modes are calculated numerically, presented, and discussed. We 
conclude in Section |5] with a speculative discussion of the impli- 
cations our results have for the long-term thermal and dynamic sta- 
bility of cool-core galaxy clusters. 



the first and second adiabatic invariants of each particle. These en- 
sure that any change in magnetic field strength B and/or density 
p must be accompanied by corr esponding changes in and py 
dChew. Goldberger & Lowlll956l) . 

The pressure anisotropy impinges noticeably on the buoyancy 
instabilities that afflict stratified plasmas, such as the magnetother- 
mal instability (MTI) and the HBI. Because the anisotropy impedes 
the convergence and divergence of magnetic field lines, the MTI 
mechanism is strengthened, while the HBI is suppressed. In partic- 
ular, pressure anisotropy shifts the fastest-growing HBI modes to 
much longer wavelengths than the very short scales it was thought 
to favour. In fact. ,Kuna t2011i) finds that the fastest growth occurs 
for a wavelength ~ 10(Amfp^^)^'^'^, which in the weakly collisional 
ICM is approximately global, Q.2H-H. These longer modes op- 
erate on a timescale slower than conduction but faster than viscous 
diffusion; thus they can effectively access free energy via the heat- 
flux, while minimising viscous damping. 

Pressure anisotropics also give rise to a host of 'microinsta- 
bilities', which include the firehose, mirror, and gyrothermal insta- 
bilities ( Schekochihin et al. 2005, 2010). These generate turbulent 
fluctuations on 'nano-scales' of ~10 npc and ~10 hr, which likely 
set the pressure a nisotropy and heat-fluxes on the larger scales that 
the HBI inhabits dSchekochihin et alj|2010l) . Ideally, these fluctua- 
tions would be 'smoothed out' by a mean field theory or else ac- 
counted for in some fashion. In this work we assume from the out- 
set that our plasma is Maxwellian, and so these instabilities do not 
appear in our linear theory (though they may during the nonlinear 
phase of the evolution). 

2.2 Governing equations 

The model we use is that of Braginskii-MHD, in which the equa- 
tions of classical MHD are employed but with speci al prescrip- 
tions (closures) for the transport of momentum and heat dBraginskiil 
[l963). The evolutionary equations governing the mass density p, 
velocity v, dimensionless entropy S — \n[p~^^^p), and magnetic 
field B are, respectively. 
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2 FORMULATION OF THE PROBLEM 
2.1 Plasma dynamics 

The macroscopic scales of the ICM (e.g. the thermal-pressure scale 
height H) are much greater than the particle mean free path Amfp 
and thus a (collisional) MHD description of the ICM plasma re- 
mains valid. On the other hand, A^fp is much larger than the ion 
gyroradius, a physical regime that ensures the predominance of cy- 
clotronic motions. As a consequence, the transport properties of 
the plasma deviate significantly from what we might expect from 
classical MHD. For instance, the conduction of heat is strongly 
anisotropic with respect to the local magnetic field direction, a 
property that fundamentally alters the ICM's convective stability 
(Balbus 2000). Temperature gradients, rather than entropy gradi- 
ents, become the discriminating quantities t hat determine stabil ity, 
regardless of whether temperature increases ( lBalbusl200dl200lh or 
decreases ( IOuataerj|2008l) in the direction of gravity. In addition, 
the fluid is subject to pressure anisotropics, i.e. the gas pressure 
perpendicular p± and parallel py to the local magnetic field may 
not be equal. Pressure anisotropics arise from the conservation of 
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is the convective (Lagrangian) derivative, g is the gravitational ac- 
celeration, and P is the (thermal) pressure tensor with the isotropic 
pressure defined through p — (l/3)TrP. The viscous heating 
rate is F, the radiative cooling rate is A, and q is the conductive 
heat-flux. In addition, the magnetic field must satisfy V-B — 0. 
Throughout we assume that the gas satisfies the ideal gas law 
p — pv'ti,, where v'^y^ = 2kB.T /nip is the thermal speed of the 
ions (with fee denoting the Boltzmann constant and rrip proton 
mass). The adiabatic index 7 has been set to 5/3. We consider a 
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hydrogenic plasma with equal ion and electron number densities, 
m = We, and temperatures, Ti = Te = T. 

2.2.1 Pressure tensor 

The ICM plasma distribution function is gyrotropic. As a conse- 
quence, the pressure tensor reduces to: 



P=px(l-bb)+pi,bb, 



(6) 



where b — B /B, and px and are the pressures perpendicular 
and parallel to the local magnetic field. The total gas pressure is 



P ■ 



2 1 

3PX + -P11 



(7) 



When the ion-ion collision time m is much shorter than the charac- 
teristic timescales associated with the macroscopic fields, the pres- 
sure anisotropy may be computed from 



p± ~P\\ = 0.960 piTii — In — - . 
" at 



(8) 



Here we have ignored the contribution of the electrons, which 
is a factor ^(m.e/m iY^'^ smaller than that of the ions 
dCatto & Simakovll2004h . By using equations lO and l|4) to replace 
the time derivatives of density and magnetic field strength with ve- 
locity gradients, equation ([Sj may be written as 
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= 3pv (bb- -\) -.Vv 
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where we have introduced the (kinematic) viscosity coefficient 
u = 0.960 X ^v^kTii, 
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This pressure anisotropy is the physical effect behind what is 
known as Braginskii ( 1965) viscosity - the restriction of the vis- 
cous damping to the motions and gradients parallel to the magnetic 
field. 

Other than impeding motions in the momentum equation l|2]l, 
the Braginskii stress appears in the entropy equation l[3) through the 
viscous heating term, F. Using the closure l[8j, this heating term is 
oc !-'ii(p± and is therefore higher order in the linear analysis 

we perform. 



2.2.2 Heat conduction 

When the particle gyroradius is much smaller than the collisional 
mean free path, heat is res tricted to flow along magnetic lines of 
force (e.g. lBraginskii|[l96^ . The heat-flux can then be written as 



(11) 



The parallel thermal diffusivity k is dominated by the contribution 
from the electrons. 
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where Wt^.c = 2k^T /rric is the thermal spee d of the electrons 
and T ee is the electron-electron collision time dCatto & Simakovl 
|2004) . The parallel thermal diffusivity of the ions is a factor 



~(me/mi)^/^ smaller jS pitzeJ 19621) . For future reference, we also 
define the (Spitzer) parallel thermal conductivity 

X = (P/T)^^. (13) 
2.2.3 Radiative cooling 

The second term on the right side of (O is the radiative cooling rate 
function A. We adopt thermal Bremsstrahlung radiation, and take 
A to be 
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jRvbicki & Lightman|[T97^ . The cooling in the ICM is dominated 
by Bremsstrahlung above temperatures ~1 keV. 

2.3 Model geometry 

The stability calculation requires us to set up a suitable geometry 
and background magnetic field configuration. Galaxy clusters are 
approximately spherical, but in this geometry it is unclear what the 
most appropriate (or analytically feasible) configuration the mag- 
netic field should take. Previous linear theory has sidestepped these 
issues by employing a local model, in which the magnetic field 
may take a simple form (purely vertical, for instance). However, 
in this paper we explicitly wish to test the importa nce of globa l 
structure when in the presence of Braginskii viscosity ( lKunzl20Tlh . 
We hence employ an intermediate model, the plane-parallel atmo- 
sphere (se e e.g. .Lamb 1997b, u s ed frequently in stellar co nvection 
problems jHurlburt et alj|l989l; [Brandenburg et ai]|l996l) . This is 
an idealised physical system, but one that can help isolate features 
associated with a cluster's global radial structure while evading the 
uncertainties and analytic difficulties of magnetic field configura- 
tions in spherical geometry. Though our model poorly describes 
a cluster in its full complexity, it should nevertheless exhibit its 
salient physics, shared by more advanced models. 

The plane-parallel atmosphere is an approximation to a chunk 
of the spherical cluster. Terms arising from spherical geometry are 
dropped; angular structure is neglected but the radial structure is 
retained. Put another way, the model is local in the angular (i.e. 
horizontal, x, y) directions and global in the radial (i.e. vertical, z). 
We adopt a uniform gravitational acceleration in the vertical direc- 
tion, g = —gz. The vertical extent of the layer is taken to be finite 
and conditions need to be applied at its upper and lower bound- 
aries, which occur at z = and z = Z. Finally, the background 
magnetic field is assumed to be constant and vertical within the 
layer. The plane-parallel atmosphere is well-suited to the morphol- 
ogy of the HBI modes in this set-up beca use these p ossess small 
radial and large horizontal wavenumbers ( ,Kunzll2oT Ih. The ' quasi- 
globar plane-parallel atmosphere works best in the outskirts of the 
cool core: neglecting the geometric spherical terms introduce er- 
rors of order Z/r. However, these errors we believe will not alter 
the qualitative behaviour revealed by our calculations. 



3 EQUILIBRIUM STATE 

We consider a static Maxwellian plasma vertically-stratified in both 
density p[z) and temperature T{z), threaded by a uniform back- 
ground magnetic field oriented along z. We further assume that the 
ratio of the thermal and magnetic pressure is large: 

Sttp ^ 2i;t\ 
52 



> 1, 



(15) 
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where va = B /{^■npY^'^ is the Alfven speed. Faraday rotation 
measurements suggest j3 ~ 10^-10^ in the co ol cores of galaxy 
clusters (for a review, see lCarilli & TavloJ2002h . 

The equilibrium must satisfy force and energy balances 



dp J 

— = —an and — — —A, 

where the heat-flux in the background state is given by 



dT 
Az • 
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We consider two thermal equilibria: one in which there is no cool- 
ing, A = 0, and hence no conductive heating; and one in which 
A 7^ and Bremsstralung radiation balances conductive heating. 

Finally, thermal boundary conditions need to be applied at the 
upper and lower boundaries of the layer, at z = and z = Z. 
Various boundary conditions may be invoked, but we choose the 
simplest: we force the top and bottom layers to have fixed temper- 
atures: T(0) = To and T{Z) = Tz- 



3.1 No cooling: A = 

In order to preserve thermal equilibrium when cooling is absent, 
the background heat-flux must be constant: 



(18) 



Enforcing the boundary conditions, equation dlSb may be inte- 
grated to yield the temperature profile 

. 2/7 



T{z) = To (l + C J 



(19) 



where = [(T^/Tq)^''^ — 1] measures the magnitude of the con- 
stant heat-flux through the atmosphere (g oc —Q. Combining this 
result with equation l ll6t determines the pressure profile 



where 

Ho = Vti,^o/g 



5/7 



- 1 



(20) 



(21) 



(22) 



is the thermal-pressure scale height at 2 = 0, and po is the pressure 
at 2 = 0. The density follows from the equation of state. 

This equilibrium is characterized by two dimensionless pa- 
rameters: ^, which is ~ 10-50 in typical cool-core clusters, and 

which is a measure of the vertical extent of the layer One example 
of this atmosphere with Tz /To = 2.5 (i.e. ( ~ 23.7) and C/ = 2 is 
shown in Figure[T]as the red lines. 

3.2 Bremsstrahlung cooling: A oc p^T^/^ 

When Bremsstrahlung cooling is present, no analytic solution ex- 
ists. Instead, we employ a shooting method to solve equations U6t- 
l ll7t , subject to three boundary conditions: T — To and p — po at 
= 0, and T — Tz al z = Z. The presence of cooling introduces 
another dimensionless free parameter. 
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Figure 1. Equilibrium atmospheres with (blue lines) and without (red lines) 
Bremsstrahlung cooling for Tz/To = 2.5 (i.e. C ^ 23.7) and g = 2. 



which is the ratio of the radiative diffusion time (across a distance 
Z) and the characteristic cooling time (each evaluated at 2 = 0). 
S is related to the Stefan number (Sf) used in some convection 
studies. If we scale space by Z, phy po, and T by To, the equations 
can be written in the following form, suitable for numerical work: 



Aq 
Az 



-Sp^T 



1/2 



dT 5/2 
Az 



^^-gpT-'+pqT-'/\ 
Az 



(24) 



(25) 



(26) 



(23) 



The boundary conditions are hence p{Q) = 1, T(0) = 1 and 
T{Z) = (Tz/To). Note that in the limit of small S (ineffective 
radiative cooling), the first equation gives a constant q to leading 
order and we recover the 'no-cooling' solution of the previous sub- 
section. 

Figure [T] describes an example of a cooling atmosphere with 
Tz/To = 2.5 (i.e. C 23.7), = 2, and <S = 45, plotted with 
the blue lines. This atmosphere i s very similar to the co ol core of 
A85, particularly for 2 > 50 kpc dCavagnolo et alj200^ . Note that 
the effects of cooling are localised to relatively small z, where the 
density is greatest. As a result, the heat conduction q (not plotted) is 
attenuated at these 2, dropping to approximately zero at 2 = 0. The 
gradient of the temperature follows, via and is hence flatter 
than in the non-cooling atmosphere. 



4 LINEAR MODES 
4.1 Linearised equations 

We consider two-dimensional perturbations, 5p, ST, SB, and Sv, 
upon our equilibrium state that exhibit a space-time dependence 
oc f{z) exp{at + ikx), where a is a (complex) growth rate and k 
a (real) wavenumber Such modes are local in the direction perpen- 
dicular to the background magnetic field (a;) and global in the direc- 
tion along the background magnetic field (z). Because the modes 
are two-dimensional, it is convenient to use the magnetic flux func- 
tion A defined through B — Vx(Ay). Finally, we scale space by 
Z, speed by Vth,o, density by po, temperature by To, and the flux 
function by BoZ. Time, as a consequence, is scaled by Z/vth,o- 
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The governing equations, written to linear order in the pertur- 
bation amphtudes, are then 



5p ... / d In p 9 , r 
a — = -\kdvx - I 92 ' 
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(30) 
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These equations introduce tliree dimensionless parameters: the 
plasma beta at z = 0, denoted by /3o, and the Reynolds and Peclet 
numbers at z = 0, defined through 



„ Zvth,o J n Zvth.o 
Re — — and Pe — 



(33) 



which quantify the relative importance of viscous and thermal 
transport, respectively. In the context of the ICM, both numbers 
can be replaced by a single parameter, the inverse Knudsen number 
at « = 0, defined through 

Ho 



Knn ' = 



Amfp.O 

1207 



10~* cms"2 



-1 (^K^\ 
Vo.Ol cm-3/ 



(34) 



2keV 



The Knudsen number is simply the ratio of the mean free path to 
the scale height. We find that 



Re = 2.08eKno^ and Pe = 0.042 ^KnoV 



(35) 



There are hence six dimensionless parameters that fully specify the 
problem. The equilibrium is set by ^, Q, and S, while the linear 
analysis introduces Po, Kno, and the (scaled) horizontal wavenum- 
ber k. 

Finally, the Prandtl number may be written as 

1/2 

~ 0.02. (36) 



Pr = 0.607 



Though Pr does not appear explicitly in the governing equations it 
is still of interest. Note that Pr is roughly constant. This implies that 
viscous forces operate on a timescale TvIsc that is a fixed number 
larger than the timescale on which conduction operates, Tcond- For 
the approximately incompressible pert urbations ex amined in this 
paper, Tvisc/rcond = (2/15)Pr-^ ^ 6 ( lKunj|201 ih . 



4.1.1 Boundary conditions 

This set of coupled ordinary differential equations must be solved 
over the domain z = [0, 1] subject to boundary conditions. We 
first assume that there is no penetration through the upper and lower 
boundaries of the atmosphere, and so Svz = at z = and 1. In 
addition, the magnetic field is constrained to be purely vertical at 
the boundaries: dSA/dz = at 2 = 0, 1. Now, if we enforce 
constant temperatures To and Tz at z — 0, 1, then we must also 
have zero temperature fluctuations at the boundary, ST = 0. This 
boundary condition is used in most numerical studies of the HBI 
in plane-parallel geometry and is the simplest to implement. Other 
conditions were trialed, such as zero heat-flux perturbation on the 
lower boundary, with little change in the qualitative results. 



4.2 Numerical scheme 

We solve equations ll27t-l|29ll numerically via a pseudospectral 
technique. The z-domain is partitioned into a Gauss-Lobatto grid 
with A'^ = 250 grid points and the differential operator is dis- 
cretized as a (Chebyshev) derivative matrix (Boyd 2000). A 5A'' x 
5A'^ matrix equation ensues taking the form of a generalised alge- 
braic eigenvalue problem; this yields a spectrum that approximates 
that of the original differential operator. The eigenvalues a are ob- 
tained using standard linear algebra routines such as the QZ algo- 
rithm or an Arnoldi method (Golub & van Loan 1996). 

4.3 Results 

We present results for various parameters. Throughout, however, 
we set g = 2, Tz /To = 2.5 (or equivalently C = 23.7). This 
leaves the cooling parameter S, which we either assign to be or 
45; the plasma beta at the lower boundary /?o = 10^; and the in- 
verse Knudsen number Kuq^, which we take usually to be 1500. 
Finally, the horizontal wavenumber k ranges between and 1000. 



4.3.1 Growth rates 

On account of the finite vertical domain of the layer, there are a 
discrete number of HBI modes each with a distinct vertical struc- 
ture, as opposed to the kz continuum exhibited by the local model. 
The fastest growing of these modes undergo the least vertical vari- 
ation (corresponding roughly to smaller k^ in the local analysis), 
while the slowest growing modes exhibit much finer-scale struc- 
ture (larger kz). 

In Fig. |2] the growth rates of the leading modes are plotted 
against horizontal wavenumber k for the non-cooling and cool- 
ing cases. In accord with the local analysis, we find that the 
fastest-growing HBI modes have large horizontal (perpendicular) 
wavenumber, with kZ > SQKtiq^^^ (see figure 3 of lKunzlbol ih . 
This wave configuration maximizes convergence/divergence of the 
background heat-flux and the consequent heating/cooling of the 
plasma. Modes on arbitrarily small horizontal scales grow at the 
maximum rate, and so the system is ill-posed unless these small 
scales are regularised via the inclusion of finite Larmor radius 
(FLR) effects, for example. Conversely, most modes are extin- 
guished for sufficiently small k (long horizontal wavelengths): an 
HBI mode ceases to grow efficiently when its horizontal wave- 
length becomes much longer than its characteristic vertical scale. 
The concentration of the heat-flux becomes inefficient in this case. 

The magnitude of the growth rates that appear in Fig.|2]is typ- 
ical for sensible parameter choices. The peak e-folding times are 



©2011 RAS, MNRAS 000,[T]{To] 



6 Latter & Kunz 




1000 



Figure 2. Growth rates of various HBI modes as functions of the perpendicular wavenumber k. The parameters of the system are: Tz /To = 2.5, Q = 2, 
Po = 10^, and Kiip ^ = 1500. The left panel con'esponds to the non-cooling case 5 = 0, the light panel to <S = 45. The growth rates are either purely real 
(blue lines) or come in complex conjugate pairs (red lines). Note that there exist multiple bifurcations at which conjugate pairs 'detach' and become purely 
monotonic modes, or two monotonic modes coalesce and form a conjugate pair The cooling case appears to favour oscillatory growth (complex conjugate 
pairs), especially for large k. 



roughly Z/i;tii,o « 400 Myr for Z = 250 kpc and fcBTo = 2 keV. 
As /3o is lowered, slower-growing modes, which exhibit finer verti- 
cal structure, are stabilised because the Alfven length reaches their 
characteristic lengthscales and magnetic tension comes into play. 
For instance, when kZ — 300 and /3o — 10^ (with no cooling) 
there are 32 growing modes; but when /3o is decreased to 10* this 
number drops to 14, and then to 5 when /3o — 10'^. Greater mag- 
netic tension (lower /3o) also lowers the maximum growth rate. 

Unstable HBI modes can exhibit either monotonic or oscilla- 
tory growth analogous to magnetoconvection in a finite layer (e.g. 
IProctor&WeisJ 19821) . In the case of oscillatory growth, modes oc- 
cur in complex conjugate pairs. As k (or Kno) varies, bifurcations 
occur where conjugate pairs 'detach' and become purely mono- 
tonic modes, or monotonic modes coalesce and transform into con- 
jugate pairs. In Fig. |2] these two modes are represented by either 
blue (purely real) or red (complex) lines. This behaviour should be 
contrasted to the local analysis in which the HBI is always mono- 
tonically growing. It follows, in part, from the variation in the back- 
ground temperature and, in particular, its influence on the thermal 
conductivity. In a global model the perturbed heat flux introduces 
a term proportional to {dT/dz)5x oc {dT/dz)5T, which is absent 
from a local analysis. Under the right conditions, this term can fa- 
cilitate growth by injecting the free energy from the gradient into a 
thermodynamic wave during its 'compressed' phase. Interestingly, 
in a radiatively cooling atmosphere growing modes tend to favour 
the oscillatory form. Cooling accentuates local perturbations in 5T 
and thus reinforces the oscillatory response of the perturbed heat- 
flux. 



4.3.2 Eigenfunctions 

In Fig. [3] we show the four fastest-growing eigenmodes for <S = 0, 
Kn^^ = 1500, and /3o = 10^. The horizontal wavenumber 
k = 250. We plot their Lagrangian displacements, and ^z, the 
scaled density Sp/p and temperature 5T/T, as well as the per- 
turbed magnetic-field lines. 

What is most striking in these plots is that the the fastest- 
growing modes, which should dominate the driving of nonlinear 



disordered flows, are confined to the innermost ~20% of the at- 
mosphere, as presaged by the local WKBJ approach (Kunz 2011). 
There are two reasons for this. First, because the HBI growth rate is 
proportional to {dlnT/dlnz)^^^, these modes take advantage of 
the steeper temperature gradient at small z. Indeed, the local growth 
rate is a factor C^^^^ ~ 5 less at the top of the layer than it is at the 
bottom. Second, because the kinematic viscosity is u oc 
the higher temperatures and lower densities characteristic of the up- 
per portion of the atmosphere guarantee rapid viscous damping of 
wavelengths significantly smaller than the local thermal-pressure 
scale height. This is also the cause of the relatively straight field 
lines for z > Q.2Z for all of the shown eigenmodes. As a con- 
sequence of these two effects, the degree of confinement depends 
on the parameters (Tz/Tq) and Kno. We find that the fraction 
of the atmosphere perturbed by the fastest growing mode scales 
like {Tz/To)~^^^Knl/^ for Kno > 1/3000. The dependence on 
Kno here coincides, unsurprisingly, with the HB I's favoured verti- 
cal wavelength in the local analysis jKunzll20Ilh . 

Slower growing modes (not shown) exhibit greater variation 
(i.e. more nodes) and extend farther into the atmosphere. Yet only 
the very slowest modes exhibit significant perturbations throughout 
the entire layer. For example, a typical e-folding time for the growth 
of substantial magnetic structures at large z is >5 Gyr, an order of 
magnitude longer than the peak e-folding time. Consequently, it 
is unlikely that such modes are dynamically important during the 
nonlinear saturation of the HBI. 

Figure|4]presents the four fastest-growing eigenmodes for the 
same parameters as in Fig. [3] except that S — 45. These modes 
occur in complex conjugate pairs, but we plot only one member 
of the pair. The modes' morphologies are qualitatively the same 
as in the non-cooling case; however the localisation is less pro- 
nounced. This is due in part because the equilibrium temperature 
gradient and heat-flux are smaller at low z when cooling is present. 
Thus both the free energy source and the catalyst for instability are 
diminished. To compensate, the mode 'spreads out' to access as 
much energy as viscous damping permits. Nonetheless, maximum 
HBI growth rates suffer and drop, in relative terms, by roughly a 
third (cf. Fig.[2j. As with the non-cooling atmosphere, the growth 
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Figure 3. Eigenvectors (^a; = &Vx/cr, (^^ = 5vz/cr, Sp/ p, 5T/T) and magnetic field lines of the four fastest-growing modes for 5 = 0, Kn^^ = 1500, 
/3o = 10^, and the horizontal wavenumber kZ = 250. The real (imaginary) part of each eigenvector is denoted by the solid (dashed) line. The growth rates 
are given at the top of each plot; for reference, Z/vh^ q a: 400 Myr for Z = 250 kpc and /cbTq = 2 keV. 



of significant magnetic structures at large z is relatively slow, with 
e-folding times of at least 5 Gyr. 



4.4 Comparison with simulations 

Other than elucidating the nature of the HBI under global clus- 
ter conditions, the above analysis also serves as a useful testbed 
for numerical codes that include anisotropic conduction and vis- 
cosity. In this subsection we offer a comparison of non-radiative 
HBI growth rates computed from the linear theory (Section [4.3.1t 
and from simulations using the Godonov code Athena ( Stone et al. 
I2OO8 ) equipped with a Braginskii stress and anisotropic heat con- 
duction. Specific ally, we use the s imulation run 'H2dBrag' pre- 
sented in § 4.3 of lKunz et all ilOlj) . 

The simulation was initialized with an atmosphere identical to 
that denoted by the red lines in Fig.[T]with parameters the same as 
in Fig.[3l i.e. /3o = 10^ and Kn^ ^ = 1500. Likewise, the bound- 



ary conditions on the top and bottom of the computational domain, 
which employed 512x1024 zones to span Ho x 2Ho, were the 
same as those described in Section l4. 1 . 1 1 The equilibrium state was 
broken by seeding the velocity with small-amplitude white noise. 
Further details reg arding this simulation may be found in §§ 3-4 of 
iKunz et al.l l2oT2h . 

In Fig. |5] we superimpose the numerical dispersion rela- 
tion (asterisks) on the analytical dispersion relation for the two 
fastest-growing modes. The numerical growth rates were obtained 
by Fourier transforming a periodic reconstruction of the simula- 
tion domain, selecting prominent and easily identifiable modes, 
and tracing their evolution throughout the linear phase of growth 
(tVth.o/Z ~ 2-5). The agreement between the code results and 
the analytical prediction is extremely good, up to wavenumbers ap- 
proaching those associated with the grid spacing (fcmax = 10247r). 
The corresponding real-space profiles of the temperature and den- 
sity fluctuations (not shown here) are also in remarkable agreement 
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Figure 4. Eigenvectors (^a; = &Vx/cf, = Sv^/cr, 5p/ p, 5T/T) and magnetic field lines of the four fastest-growing modes for a cooling atmosphere. 
Parameters are the same as in Fig.|3]except that S = 45. These modes occur in complex conjugate pairs, but we plot only one member of each pair 



with those predicted in Fig. |3] This is a valuable check and cross- 
validates both the linear theory and the numerical code. 

Because of numerical dissipation on the grid, at fc > 
the numerical growth rates abruptly depart from the theoretical pre- 
diction and drop precipitously. At which k this 'unphysical' regu- 
larisation occurs depends entirely on the resolution adopted, and 
formally such simulations can never be resolved. In practice, how- 
ever, it is unclear how many of these small-scale HBI modes are 
actually required to adequately describe the nonlinear dynamics of 
the system. Also unclear is whether the grid regularisation realisti- 
cally approximates regularisation by FLR effects, and if this mat- 
ters in the saturat ed state. Some of these issues are addressed in 
iKunzetalldioT^ . 



5 DISCUSSION 

In this paper we have examined the linear stability of weakly colli- 
sional, thermally stratified atmospheres in which the conduction of 



heat and the viscous dissipation of motions are anisotropic with re- 
spect to the magnetic field. Such atmospheres are representative of 
those found in the cooling cores of non-isothermal galaxy clusters, 
which are subject to a heat-flux buoyancy-driven instability (HBI). 

Numerical simulations of the HBI have established that the 
instability acts in such a way as to disrupt the conductive flow of 
heat, ultimately insulating the core from the hot thermal bath at 
large radii. In the presence of appreciable radiative losses, a cooling 
flow inevitably develops unless turbulent stirring by other means 
can reopen the field lines and reinstate the conductive flux. 

The pressure anisotropies generated by the HBI play an im- 
portant role in suppressing the instability by reducing growth rates 
and by shifting the fastest-growing modes to relatively long wave- 
lengths. In the bulk of cool-core clusters, the collisionality is suf- 
ficiently low that these wavelengths become comparable to the 
thermal-pressure scale height of the gas. The HBI then operates 
on a global scale, and the analysis we have performed here marks 
one contribution towards understanding how the HBI changes in a 
global setting. 



© 2011 RAS, MNRAS QOQ.mflOl 



The HBI in a quasi-global model of the ICM 9 









1.4 








- 


/ % ; 








1.2 


- 










1.0 


- 


// 


0.8 


- 

- 


/ 

1/ - 


0.6 










* run H2dBrag of 






Kunz et al. (2012) ' 


0.4 




1 1 



10 100 1000 

kZ 

Figure 5. Comparison between the analytical dispersion relation (red and 
blue lines; as in Fig. [3) and t he numerical dispe rsion relation (asterisks; ob- 
tained from run HldBrag of jKunz et al.ll2012l) for the two fastest-growing 
modes. The parameters of the system are: T^/Tq = 2.5, Q = 2, 
/3o = 10^, and Kn^ ^ = 1500. 



Our main result is that the low coUisionality beyond ~50 kpc 
in typical cool-core clusters significantly impedes the otherwise 
disordered motions driven by the HBI. As a result, the magnetic- 
field lines can remain relatively straight over a significant fraction 
of the Hubble time and the conductive flow of heat may proceed at 
an appreciable fraction of the Spitzer value. Near-complete field- 
line insulation due to the HBI could very well occur, but our anal- 
ysis suggests that it is well-localised to the innermost regions of 
the cool core. This highlights the need for radio-mode feedback at 
these scales from a powerful central dominant galaxy. 

This conjecture is supported by rece nt numerical simulations 
of the HBI presented in lKunzetalj32012h . which self-consistently 
take into account the pressure anisotropy driven by the HBI and 
its a dverse consequences for the instability's development. In con- 
trast, [p^risheTd] ([lOl^, hereafter Pi 2) perform analogous simu- 
lations of the HBI but these show that pressure anisotropy has little 
effect on the nonlinear saturation of the instability; the HBI op- 
erates throughout the bulk of the cluster core unhampered and im- 
poses significant field-line insulation throughout. These very differ- 
ent conclusions possibly arise from dissimilar parameters and mod- 
els. In particular, the PI2 simulations that employ a plane-parallel 
atmosphere are of very small vertical extent, typically Z = 0.2Ho- 
As a consequence, there is insufficient room for global effects to 
come into play. These simulations cannot exhibit the mode locali- 
sation that we emphasise. However, P12 also perform two spherical 
global simulations of the HBI in model cool-core clusters, and these 
support the results of their plane-parallel simulations. Here the con- 
tradiction with our work is more puzzling, but probably arises from 
different parameters. For instance, P12's ICM model is far more 
collisional and less magnetised than ours, with /3 ~ 10''. The ef- 
fective absence of magnetic tension in these simulations may be 
important as magnetic tension suppresses those shorter-scale (lo- 
cal) HBI modes that function throughout the core. Alternatively, 
the discrepancy may be due to P12's different initial magnetic-field 
geometry (tangled on scales 30-50 kpc) or due to a basic dis- 
similarity between the spherical and plane-parallel geometries. But 
the apparent disagreement is a spur for future work: further three- 
dimensional simulations of cool-core clusters under representative 
conditions are needed. 
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